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Disordered biopolymer gels have striking mechanical properties including strong nonlinearities. 

In the case of athermal gels (such as collagen-I) the nonlinearity has long been associated with a 
crossover from a bending dominated to a stretching dominated regime of elasticity. The physics 
of this crossover is related to the existence of a central-force isostatic point and to the fact that 
for most gels the bending modulus is small. This crossover induces scaling behavior for the elastic 
moduli. In particular, for linear elasticity such a scaling law has been demonstrated [Broedersz et 
al. Nature Physics, 2011 7, 983]. In this work we generalize the scaling to the nonlinear regime with 
a two-parameter scaling law involving three critical exponents. We test the scaling law numerically 
for two disordered lattice models, and find a good scaling collapse for the shear modulus in both 
the linear and nonlinear regimes. We compute all the critical exponents for the two lattice models 
and discuss the applicability of our results to real systems. 


I. INTRODUCTION 

Collagen is the single most abundant protein in the 
animal kingdom; most of it occurs in the form of fibrous 
collagen such as collagen-I. A major component of tissue 
is disordered networks of fibrous collagen, i.e. collagen 
polymer gels. These networks are one of many impor¬ 
tant biological examples of “athermal” biopolymer gels - 
athermal means that the coiling of the stiff fibrils due to 
thermal fluctuations is negligible. These materials have 
remarkable mechanical properties: even though the con¬ 
stituents are in the linear elastic regime (up to the yield 
point of the network), the network itself is nonlinear and 
strain-stiffening. For example, in dilute collagen-I gels 
the shear modulus increases by orders of magnitude as 
the strain increases [1-5]. Such nonlinear effects are im¬ 
portant in biology, for example in the stiffening and align¬ 
ment of tissue near a growing tumor [6]. 

The unusual mechanical properties of collagen are be¬ 
lieved to be intimately connected to two facts. First, the 
hbrils have a large stretching modulus, k, and a small 
bending modulus, k. More precisely, we define R = n/ka^ 
where a is the mesh size of the network; typically k <C 1. 
For example, for collagen-I we can think of the fibrils 
as elastic rods. In the dilute case (densities of order 
2mg/ml) R ^ {d/a)^ ~ 10“^ where d ~ 30nm is the 
diameter of the fibrils and a ^ 1^ is the mesh size. Sec¬ 
ond, the disordered structure of collagen can accommo¬ 
date small deformations through bending the fibrils, but 
at large deformation, the bending modes are exhausted 
and the fibrils have to be stretched. Owing to the huge 
difference between the bending and stretching stiffness of 
the fibrils, such a bending-to-stretching crossover lead to 
a dramatic increase of the elastic moduli [7-11]. 

The existence of such a bending-to-stretching crossover 
is controlled by the physics of the central-force isostatic 
point (CFIP). The concept of CFIP traces back to J. C. 
Maxwell [12]. He pointed out that the onset of mechani¬ 
cal instability for a system of particles with central force 


interactions occurs when the mean coordination num¬ 
ber of bonds is twice the spatial dimension, (z) = 2d. 
At this point the number of degrees of freedom, d, and 
the number of constraints, 2:/2, of each particle are the 
same [13-20]. For (z) < 2d there are soft floppy modes; 
for (z) > 2d there are redundant bonds and states of self¬ 
stress. Thus the CFIP is a critical point of a mechanical 
phase transition. Biopolymer gels with k <g; 1 can be 
thought of as a nearly central-force system if we identify 
the crosslinking points as the particles and the polymer 
fibrils between them as bonds. The coordination number 
z depends on the type of gel. Each crosslinking point can 
have a maximum coordination number Zmax, depending 
on the type of bonding. The most common situation, two 
polymers meeting at a crosslinking point, corresponds to 
-Zinax = 4 (Collagen I belongs to this class). The aver¬ 
age coordination number (z) < Zmax because of dangling 
ends. Therefore, pure Collagen I and similar biopolymer 
gels have mean coordination (z) < 2d and are always 
below the CFIP. In contrast, if a gel has crosslinking 
points with Zmax > 2(i the gel can go above the CFIP 
when the filaments are long. This may be the case for 
multi-component complex gels seen in extracellular ma¬ 
trix material. 

Since the CFIP is a critical point, it is natural to sus¬ 
pect that mechanical properties such as the shear mod¬ 
ulus would exhibit scaling properties in its vicinity, and, 
indeed, for linear elasticity such scaling has been demon¬ 
strated [4, 21, 22]. In this paper we extend scaling re¬ 
lations near the CFIP to the nonlinear regime. As sim¬ 
plified models of real biopolymers we use simulations on 
two different disordered lattice models, the diluted tri¬ 
angular and kagome lattices (see Fig. I (a,b))[4, 21, 22]. 
These two lattices have maximum coordination Zmax = 6 
and Zmax = 4; they represent two classes of gels with 
^^max > 2d and Zmax < 2d. We use these lattices to verify 
our scaling laws and to extract exponents. We discuss 
the application of the scaling regimes to more realistic 
biopolymer networks in Sec. IV. 
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FIG. 1. Examples of network configurations before (grey) 
and after (red) a large deformation for the diluted triangular 
lattice at k = 10“® (a) and kagome lattice at k = 10~^. 
(b). Shear modulus as a function of strain for varying values 
of p (see legends) for the diluted triangular lattice (c) and 
kagome lattice (d). Point A(B) corresponds to the gray(red) 
configuration in Fig lab. 


The scaling form that has been proposed [4, 21, 22] 
previously for the linear shear modulus, G is: 

GiAz, k) = k\Az\f e±, linear {fi/\Az\^) ■ (1-1) 


Here Az = (z) — 2d is the extra coordination above the 
CFIP and is a scaling function. The exponent / cap¬ 
tures the emergence of rigidity as coordination increases 
in the central-force system, and the crossover exponent 
(j) captures the effect of k as a relevant perturbation near 
isostaticity. This scaling law for the linear elasticity in¬ 
cludes the observed effect that the elastic modulus of 
the gel can change dramatically when the connectivity is 
changed near the CFIP. Below the CFIP the elasticity of 
the gel is bending dominant, G ^ k, because it is floppy 
in the central-force limit. Above the CFIP the gel be¬ 
comes stretching dominant G ^ k. If Zmax > 2d there 
is a special regime where bending and stretching modes 
are strongly coupled, G ~ This law has 

been shown to be valid using two and three dimensional 
disordered lattice models [4, 21, 22]. 

We now turn to the nonlinear elasticity. A generalized 
scaling law should describe the very large strain-stiffening 
that is characteristic of these materials. As we discussed 
above, the physical origin of this effect is that gels in bi¬ 
ological systems are often below the CFIP {{z) < 2d). 
In the linear elasticity regime the shear modulus is pro¬ 
portional to the bending stiffness, k, and is very small 
because the shear is accommodated in the floppy bend¬ 
ing modes. As the deformation progress, the soft modes 
are stretched out and central-force paths form to bear 
the stress, increasing G\J]. 

Based on these ideas, we propose a two-parameter scal¬ 
ing law for the differential shear modulus: 

G(A0,k,7) = fclAzl^g± , (1.2) 

Here 7 is the shear strain, and the ± subscript labels the 
branches of the scaling function for Az > 0 and Az < 
0, respectively. This scaling law captures the stiffening 
effect of both the bending rigidity, k, and of the strain, 

7- 

A one-parameter version of this scaling relation for the 
nonlinear case has been previously proposed by Wyart et 
al [8]. These authors considered varying 7 with fixed k. 
They give numerical evidence that /3 = 1 for jammed 
packings, and demonstrated the scaling of the crossover 
from linear to nonlinear behavior. Our goal is to advance 
such scaling relations to include both the connectivity z 
and the dimensionless bending stiffness R as variables. 
We will collapse the whole elastic modulus curve with 
varying z and k, and analyze different elastic regimes. 
Our results are summarized in Figs. 4 and 5. 

II. MODELS AND SIMULATIONS 

The model lattices are the diluted triangular [4, 23, 24] 
and kagome lattices [15, 22], as shown in Figure lab. 
In both lattices each bond is present with probability 
p. The undiluted triangular and kagome lattices have 
coordination numbers Zmax = 6 and ^max = 4 respec¬ 
tively. It is straightforward to use Maxwell’s rule to find 
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that the CFIP is at pc = 2/3, pc = 1, respectively, at 
the mean-field level. In the triangular lattice there is a 
second-order rigidity percolation transition and the crit¬ 
ical point is slightly lower than 2/3 due to critical fluc¬ 
tuations [4, 25]. In the kagome case rigidity percolation 
is a first-order transition[22] at p = 1. For these lattices 
{z) = Zmax • (p — Pc)- These models are representative 
members of two different classes of disordered networks 
as they approach the CFIP. As above, we make contact 
with the mechanics of physical biopolymer gels by iden¬ 
tifying the lattice sites as crosslinking points and bonds 
as fibrils between crosslinks. 

We write the elastic Hamiltonian for our models as [4]: 

H = Y,\mh3f + Y. ■ ( 2 . 1 ) 

The sum (i,j) runs over all bonds and {i,j,k) runs over 
pairs of bonds that are co-linear and share lattice site j. 
We assign harmonic springs of spring constant k for the 
bonds and bending stiffness n to the bond pairs. 

In our simulations we minimize H using conjugate 
gradient methods to get the elastic energy, E. Then 
we numerically measured the differential shear modulus, 
G = at various values of p and k. Here 

A is the total lattice area. To shear the networks, the 
horizontal boundaries to which the bonds are attached 
are translated with strain 7. We use periodic boundary 
conditions on the vertical boundaries and fixed bound¬ 
ary conditions on the horizontal boundaries. Our net¬ 
work size is VF = 100 x 100. Examples of shear modulus 
data measured from our simulations before any scaling 
are shown in Fig. led. It is clear from these curves that 
strain-stiffening is prominent when the network is diluted 
(smaller p). 

In order to test the scaling relation in Eq. (1.2), we 
need to first determine the location of the CFIP and the 
exponents / and (j). For the triangular lattice, we follow 
the method of Broedersz et 0/(4] to determine the loca¬ 
tion of the rigidity percolation point for our system size, 
W. This involves looking for a peak of the nonaffinity 
parameter: 

r = ^(u. - (2.2) 

i 

in the central-force lattice (k = 0). Here i labels the sites, 
Ui is the displacement vector of the site, and is the 

displacement vector if the deformation would have been 
affine. We find that Pc = 0.651 ± 0.05 for the triangular 
lattice. Next we determine the scaling exponents / and </. 
These can be extracted from the small 7 (linear elasticity) 
limit. As shown in Figure 2 (a) we find that / = 1.3 ± 
0.1,^ = 3.6 ± 0.3 for the triangular lattice, consistent 
with previous literature on the linear elasticity of this 
lattice [4]. 

The CFIP for the kagome lattice is at p = 1, and we can 
only observe p < pc. Also, for this case, it is known [22] 




(b) 

FIG. 2. Scaling plot for linear shear modulus G (with k = 1) 
at fixed K. (a) The triangular lattice at k = 10~®. (b) The 
kagome lattice at k = 10““^. The inset shows the determi¬ 
nation of the critical point via the peak of the non-affinity 
parameter. 

that the transition is first order with pc = I; thus we 
expect, and indeed find / = 0. Our simulations show 
that <j) = 2.3 ± 0.2, as shown in Fig. 2b and consistent 
with previous literature on the linear elasticity of this 
lattice [22]. 

III. RESULTS 

To test the scaling law of Eq. (1.2) we can plot G using 
scaled parameters: 

X = k/|Az|'^, y = y/IAzj^. (3.1) 

The scaling law predicts that G/\/S.z\^ collapses into a 
single surface as a function of x, y for the correct crit¬ 
ical exponents /, </, /3. However, to verify such three- 
dimensional scaling collapse in general is computationally 
expensive. Instead, we make two types of plots using our 
simulation data. The first type of plot is at fixed k, given 
in Figure 3. We obtain the expected complete collapse 
in the nonlinear (stretching) regime where the curves are 
predicted to be independent of x. 

The second type of plot is to observe a slice of 
the scaling collapse in the three-dimensional space of 
x,y,G/\Az\^. To do this, we take a given finite value 
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FIG. 3. Scaling plot for the shear modulus at fixed R. (a) Tri¬ 
angular lattice. K — 10“®; Fig 2(b) Kagome lattice, k = 10“'*. 
As expected, collapse occurs only in the nonlinear regime. 


of a: = k/IAzI*^. In practice, this correspond to taking 
the value of R which changes as p changes so that x is 
fixed in the simulation. Recall that </> has already been 
determined from the linear regime. As a result the shear 
modulus data completely collapses for both lattices ac¬ 
cording to our two-parameter scaling, as shown in Figure 
4. In particular, for the triangular lattice G(p) lies on 
two branches, corresponding to p > Pc and p < Pc- From 
these plots we find that /? = 1.3 ± 0.1 for the triangular 
lattice and /3 = 1.7 ± 0.1 for the kagome lattice. 

We should take note of a detail in the fitting procedure 
that we and others have used [26, 27]. We are, in fact, 
never extremely close to the critical point, because suf¬ 
ficiently near to Pc there is a non-scaling feature in the 
shear modulus, a dip. We found this for both our lattice 
models, and it has been observed by others[26, 27]. It oc¬ 
curs in models which have explicit buckling[27] for fibrils 
under compression; of course, real biopolymers buckle as 
well. 

The reason is easy to see. Euler buckling is a pitchfork 
bifurcation as the compressive strain on an elastic beam 
varies. Thus we expect that the elastic energy near the 
buckling threshold, 7* will have the classic form E = 
Eo — G(7 — 7t)^i7 > 7t- The second derivative of the 
energy is the shear modulus, so we expect G —>■ G — G as 
we pass through the threshold. The background modulus 


FIG. 4. Two parameter scaling law Eq. (1.2) yields good 
collapse in both linear and nonlinear regime when x is fixed, 
(a) Triangular lattice, a;=100. (b) Kagome lattice *=0.1 


increases, so we should see a dip in G(7). 

The common way [26, 27] to introduce buckling in 
model fiber networks is to introduce an extra node in the 
middle of each link which can bend with a small bending 
modulus. There is a pitchfork bifurcation at a thresh¬ 
old just as in the case of a continuous elastic beam. In 
our models, although we do not have such buckling at 
the bond level, the crosslinking points represent internal 
degrees of freedom of the filaments and also exhibit a 
buckling instability. Nevertheless, why these show up as 
a sharp dip near pc remains a puzzle to us. 


IV. CONCLUSIONS AND DISCUSSION 

With these results in hand, we can speculate about 
their applicability to other lattice models and to phys¬ 
ical systems. Based on our results, we can discern at 
least two different classes of scaling behavior with differ¬ 
ent exponents, though both obey Eq. (1.2). We show the 
schematic phase diagrams for these two classes in Eig. 5. 

In the first case, exemplified by the diluted triangular 
lattice, we have a network that can vary through the 
CEIP as a parameter is changed. Biopolymer gels with 
-Zmax > 2(i belong to this class. Some dense forms of 
collagen, such as tendons and cartilage may also belong 
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FIG. 5. Three dimensional schematic phase diagrams for 
the triangular(a) and the kagome (b) lattices, shown in the 
space of dilution parameter p, filament bending stiffness k, 
and strain 7 . These phase diagrams show that strong strain- 
stiffening occurs when the system is bending-dominant at 
small deformation. 


to this class. In this case Zmax > 2d so that we can tune 
through Pc < 1. Then we have various regimes implied 
by Eq. (1.2): 

• A linear bending regime when Az<0,a;<Cl,?/<C 
1. This regime is described by linear elasticity so 
Q- is independent of y. This regime is bending 
dominated: 

Q-{x,y) ^ X G ^ k|Az|^“'^. (4.1) 

• A linear coupled regime when Az<0,a;^l,i/<?;1 
while K <C 1. This regime is also described by linear 
elasticity; Q- is independent of y. Also, it is in a 
critical regime (a; 3> 1) so G should be independent 


of Az. Thus: 

g- (x, y) ~ G ~ (4.2) 

• A nonlinear regime when y S> 1 while keeping 7 ^ 
1. This regime characterizes nonlinear elasticity 
as the strain goes beyond a “turning point” 7 * ^ 
I Az|^. The turning point approaches 0 as to p —Pc 
meaning that near the CFIP the linear elasticity 
regime vanishes. It is reasonable to hypothesize 
that in this regime G is independent of Az and x 
so that: 

g±{x,y)^y^^^ G ~ (4.3) 

Our simulation data verifies this hypothesis. 

• The stretching regime. This includes a small-strain 
stretching-dominated regime when Az > 0 and 
K <§; 1 , and a large-strain stretching-dominant 
regime when 7^1 and the bending modes are 
all stretched out. Here the elasticity is controlled 
by the stretching stiffness k: 

g±{x,y)I G ^ fc|Az|^. (4.4) 

The other category, for which the diluted kagome lat¬ 
tice is an example, has Zmax = 2 d and a first-order 
central-force rigidity transition at pc = 1 and / = 0 . 
In this case we have two regimes: 

• A linear bending regime when x -C I,y 1 (in 
this case we always have Az < 0) so that: 

g-{x,y) ^ X G ~ k|Az|“'^. (4.5) 

• The stretching regime, including small-strain part 
X ^ 1 and large-strain part 7 ^ 1 . In this regime 
deformations are dominated by stretching of fibers 
and 

G±{x,y)^l Gr^k. (4.6) 

The linear coupled regime and the nonlinear regime are 
missing because they are consequences of a continuous 
transition at pc < 1. As a result the crossover to stretch¬ 
ing as strain increases is more abrupt in these lattices. 

Many three-dimensional fiber-networks including di¬ 
lute collagen I, have Zmax < 2d. For this case, the central- 
force network is not rigid even at p = 1. As in the kagome 
case, scaling law predicts a linear bending regime cross¬ 
ing over to stretching at either large n or large 7 . This 
is consistent with numerical results on three-dimensional 
models [4, 28-30]. 

We can relate experimental results to our models by 
focusing on three dimensionless numbers. The first one 
is the maximum coordination of the crosslinkers z^ax 
which determines which type of the phase diagrams ap¬ 
plies. The second one is R which can be deduced, for 
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example, from the Young’s modulus and the shape of in¬ 
dividual fibers. The simplest case of this estimate has 
been given above. The third one is p. For a network 
made up of fibers of mean length L it is easy to see that 
the relation between mesh size a and fiber length L is: 
a = L/{1 — p). This allows us to estimate the effec¬ 
tive bond dilution from data. With these three numbers 
one can sketch the phase diagram and determine which 
regime the experimental system belongs to. 

In addition, it is also interesting to compare our the¬ 
ory to other types of gels with less stiff filaments and 
thus stronger thermal fluctuations. A different mecha¬ 
nism of strain-stiffening has been proposed in the litera¬ 
ture attributing strain stiffening to stretching out ther¬ 
mal undulations of polymers [3, 31]. This is relevant for 
gels composed of filaments such as actin or elastin which 
have persistence length comparable or smaller than the 
mesh size, so considerable thermal undulations present, 
thus providing a regime of entropic elasticity at small 
strain. Although for these gels stretching out thermal 
undulations contribute to the strain-stiffening, if the gel 
is below the CFIP (which is true for most bio-polymer 
gels), the bending modes still need to be stretched out in 
addition to the thermal undulations before the gel enters 


true stretching-dominant elasticity regime. 

During the preparation of our manuscript, we learned 
about two recent preprints [32, 33] which use a combina¬ 
tion of experimental and simulation tools to investigate 
nonlinear elasticity. Both of these preprints mainly fo¬ 
cus on the nonlinear shear modulus as bending stiffness 
K varies while keeping the connectivity {z) fixed. In con¬ 
trast, our study obtained full collapse of shear modulus 
in the more general case of varying both k and (z). More¬ 
over, we also identify various regimes of nonlinear elastic¬ 
ity in a three-dimensional phase diagram of connectivity, 
bending stiffness, and strain. 
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